function [rho,up,lo,bigY] = sacf(y,m,gflag)
% PURPOSE: find sample autocorrelation coefficients 
%---------------------------------------------------
% USAGE: p = sacf(y,m,gflag)
% where: y = a time-series (need not have mean zero)
%        m = # of sample autocorrelations to compute
%    gflag = 0, flag for graphing, 1 = no graph
%            (default = 0, a graph is produced)
%---------------------------------------------------
% RETURNS:
%        p = an m x 1 vector of sample acf's
%        also plots sample acf's with 2-sigma intervals
%        if gflag == 0
% --------------------------------------------------
% SEE ALSO: spacf(y,m)
%---------------------------------------------------
% Originally from teh JPL toolbox
% Vectorized by Kevin Sheppard
% kksheppard@ucsd.edu
if gflag ==0
    flag=0;
else; 
    flag=1; 
end; 
n = length(y);
rho = zeros(m,1);
npm = n+m;
tmp = std(y);
vary = tmp*tmp;
% put y in deviations from mean form
ym = mean(y);
e = zeros(n,1);
e(1:n,1) = y - ones(n,1)*ym;
bigY=repmat([1:n]',2*(m+1),1);
bigY=bigY(1:(2*n-1)*(m+1));
bigY=reshape(bigY,2*n-1,m+1);
bigY=bigY(1:n,2:m+1);
bigY=y(bigY);
E=repmat(e,1,m);
rho=(sum(E.*bigY)/(n*vary))';
up = 2*(1/sqrt(n)*ones(m,1));
lo = -2*(1/sqrt(n)*ones(m,1));
if flag == 1;
tt=1:m;
set(0,'DefaultAxesXTickLabel','remove')
set(0,'DefaultAxesXlim','remove'); 
set(0,'DefaultAxesYlim','remove'); 
set(0,'DefaultAxesXtick','remove'); 
area(rho);
title('Sample autocorrelation coefficients');
set(gca,'XLim',[1 m]); 
xlabel('k-values');
ylabel('Auto-coorelation values');
	hold on;
	plot(tt,up,'*r',tt,lo,'*r');
	hold off;
end; 